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Abstract 

Mixed linear models are commonly used in repeated measures studies. They account 
for the dependence amongst observations obtained from the same experimental unit. 
Oftentimes, the number of observations is small, and it is thus important to use 
inference strategies that incorporate small sample corrections. In this paper, we 
develop modified versions of the likelihood ratio test for fixed effects inference in 
mixed linear models. In particular, we derive a Bartlett correction to such a test 
and also to a test obtained from a modified profile likelihood function. Our results 
generalize those in Zucker et al. (Journal of the Royal Statistical Society B, 2000, 62, 
827-838) by allowing the parameter of interest to be vector- valued. Additionally, our 
Bartlett corrections allow for random effects nonlinear covariance matrix structure. 
We report numerical evidence which shows that the proposed tests display superior 
finite sample behavior relative to the standard likelihood ratio test. An application 
is also presented and discussed. 

Key words: Bartlett correction, Fixed effects, Likelihood ratio test, Mixed linear 
models, Modified profile likelihood function. 



1 Introduction 



In recent years repeated measures data have been widely analyzed in a number 
of fields, including biology and medicine. In such studies, the data are obtained 
from a number of different experimental units, each unit being observed more 
than once (Brown and Prescott, 2006). In particular, some of these studies 
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use longitudinal data (see Verbeke and Molenberghs, 2000, and Diggle et al., 
2002), in which the data are collected over time. Mixed linear models have been 
extensively used by practitioners to model repeated measures data since they 
allow the modeling of within units correlation. It is also noteworthy that there 
is available software specifically designed for the estimation of such models; 
see Pinheiro and Bates (2000) and Littel et al. (2006). For further details on 
mixed linear models the reader is referred to Demidenko (2004). 

A common shortcoming lies in the fact that in many studies the number of 
observations is small which renders approximate inferential procedures unreli- 
able. In particular, the likelihood ratio test, which is commonly used to make 
inference on the fixed effects parameters, quite often displays large size distor- 
tions when the sample size is small since its null distribution is poorly approxi- 
mated by the limiting x 2 distribution, from which critical values are obtained. 
It is possible to obtain a Bartlett correction factor and use it to modify the 
likelihood ratio test statistic in such a way to bring its null distribution closer 
to its limiting counterpart; the approximation error is reduced from 0(n _1 ) to 
0(n~ 2 ), where n is the sample size, thus making any size distortion vanish at 
a faster rate. Bartlett corrections for linear and non-linear regression models 
were obtained by Barroso and Cordeiro (2005), Cordeiro and Paula (1989), 
Cribari-Neto and Ferrari (1995), and Lemonte et al. (2008), among others. 
For a review of the literature, see Cribari-Neto and Cordeiro (1996). 

Another shortcoming relates to the effect of the nuisance parameters on the 
resulting inference on the parameters of interest. Different modifications to the 
profile likelihood function have been proposed with the aim of reducing such 
effect. The adjustment proposed by Cox and Reid (1987) can be used whe- 
never the nuisance and interest parameters are orthogonal. The Monte Carlo 
evidence reported by Ferrari et al. (2007) shows that the Cox-Reid adjustment 
can yield superior inference in finite samples. DiCiccio and Stern (1994) have 
shown that the Cox-Reid test statistic can be Bartlett-corrected, just as the 
likelihood ratio test statistic. They have extended the main result in Mukerjee 
and Chandra (1991), which had shown the Bartlett-correctability when both 
the interest and nuisance parameters are one-dimensional. The combined use 
of modified profile likelihoods and Bartlett correction can deliver accurate and 
reliable inference in small samples, as evidenced by the results in Ferrari et 
al. (2004), Ferrari et al. (2005) and Cysneiros and Ferrari (2006). 

Zucker et al. (2000) obtained improved likelihood ratio testing inference by 
deriving Bartlett corrections to the profile and modified (Cox-Reid) profile 
likelihood ratio tests for inference on the fixed effects parameters in mixed 
linear models. See also the numerical evidence in Fouladi and Shieh (2004) 
and Manor and Zucker (2004). Their results, however, are only applicable for 
testing one parameter at a time, since they only allow for a scalar parameter 
of interest. In many studies, nonetheless, practitioners wish to perform joint 
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testing inference on a set of parameters, especially when comparing three or 
more treatments in medical trials. Also, they derived the Bartlett correction 
to the profile likelihood ratio test only for the situation where the covariance 
matrix for the random effects has a linear structure. Hence, their results are not 
fully applicable in a number of situations of interest, e.g. when the responses 
of a single subject are measured sequentially and the errors are assumed to 
be correlated through an auto-correlated error structure. Our chief goal is to 
generalize their results so that they are valid in situations where the parameter 
of interest is vector-valued and the covariance matrix for the random effects 
is allowed to have a non-linear structure. We obtain the Cox-Reid profile 
likelihood adjustment, and also Bartlett correction factors for the profile and 
adjusted profile likelihood ratio test statistics. 

The paper unfolds as follows. Section 2 introduces the mixed linear model, 
Section 3 contains the three improved tests (Cox-Reid and Bartlett-corrected 
tests), and Section 4 presents numerical evidence on the finite sample behavior 
of the standard likelihood ratio test and its modified counterparts. An applica- 
tion that uses real data is presented and discussed in Section 5. Finally, Section 
6 concludes the paper. Technical details are collected in two appendices. 



2 Mixed linear models 

The mixed linear model is given by 



where y« = (yn, y i2 , • • • , yin) T is a Tj x 1 vector of responses on the ith experi- 
mental unit, (3 is an n-vector of fixed effects parameters, Xi is a r« x n known 
matrix, bj is a random effects vector (g x 1), Z, is a known Tj x q matrix, 
and 6i = (en, e^, • • • , £i Ti ) T is a T i x 1 vector of random errors. It is oftentimes 
assumed that e« ~ A/" Ti (0, cr 2 I n ), where I n denotes the Tj x Tj identity ma- 
trix and is a vector of zeros. It is also assumed that bj ~ Af q (0,G), where 
bi, b 2 , . . . , bjv, ei, €2, ■ ■ ■ ,€n are independent and G = G(g) is a q x q positive 
definite matrix, g being an m x 1 vector of unknown parameters. Model (1) 
can be written in matrix form as 



where Y = (y^, y 2 T , . . . , y^) T is T x 1, with T = ££i n , X = (Xj, Xj , . . . , 



X^) T is a T x n matrix, Z is a T x Nq diagonal matrix given by Z = 
diag(Zi, Z 2 , . . . , Zn), b = (b^, bj, . . . , b^) T is an iVg-vector and e = (ej , ej , 
. . . , eJf) T is T x 1. Thus, b ~ A/jvq(0, In<S>G), where <g> denotes the Kronecker 
product and e ~ Ar(0, o" 2 ^t); b and e are independent. 



y, = X,(3 + Z l b l + e l , i = l,...,N, 



(1) 



Y = X(3 + Zh + e, 



(2) 
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It is possible to write model (2) as 



Y = X(3 + e, (3) 

where e = Zh + e. Hence, e ~ Af T (0, E), where S = S(cj) = Z(I N ® G)Z T + 
<t 2 /t) ^ = (f? T )°" 2 ) T being an (m + 1) x 1 vector of unknown parameters. 
Hence, the log-likelihood function for model (3) can be expressed as 

£((3, W ; Y) = -| log(27r) - i log |E| - ^(Y - X/jfE" 1 ^ - X(3), (4) 
where | • | denotes matrix determinant. 

Let 9 = (tp T , <; T , uj t ) t be the (n + m + l)-vector of parameters, where ■?/> = 
(Pi, P2, ■ ■ ■ , Pp) T is the p-vetor (p < n) containg the first p elements of (5 and 
(<; T , uj t ) t is the (n — p + m + 1) x 1 vector of nuisance parameters with q = 
((3p+i, [3 p +2, ■ ■ ■ , fin) 1 ■ I n what follows we shall focus on fixed effects inference. 
In particular, we wish to test Ho : tp = against Hi : ip ^ where ip^ 
is a given p-vector. 

We follow Zucker et al. (2000) and use a reparameterization in which the 
nuisance ((<? T , ^ T ) T ) and interest (ip) parameters are orthogonal. In particular, 
we transform 9 = (tp T , ^ T , u T ) T into d = (ip T , £ T , uj t ) t , with 

£ = q + (X~ r Y,- 1 X)- 1 X~ T Y,- 1 X p iP, (5) 

where X p denotes the matrix formed out of the first p columns of X and 
X contains the remaining (n — p) columns of X. It is easy to show that ip is 
orthogonal to <fi = (£ T , u T ) T , i.e., the expected values of d 2 £($; Y)/dipd^ T and 
d 2 £(i}; Y)/dipdu)j, for j = 1,2,..., m+1, are matrices of zeros. By partitioning 
X as (X p , X) and (3 as (ip T , <? T ) T , we can write = X p ip + X^. Using (5) 
we obtain 

jcp = X> + X£, 

where X' p = [I T -X(X T 'E- 1 Xy 1 X T 'E' 1 }X p . It follows that the log-likelihood 
function in (4) can be written as 

i = eft Y) = - T - log(27r) - l - log |E| - Vs-'z, (6) 

where z = Y - X'A) - X£. 



3 Improved likelihood ratio tests 

The profile likelihood function, which only involves the vector of parameters 
of interest, is defined as £ p (ip) = £(ip, 4>(ip)), where <f>(ip) is the maximum 
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likelihood estimator of for a fixed value of if>. The likelihood ratio statistic 
for the test of Hq is 



LR = LRtyW) = 2 - W i0) )} , 

where if> denotes the maximum likelihood estimator of if}. Under the standard 
regularity conditions and under H,q, LR converges in distribution to xt- This 
first order approximation may not work well in small samples, however. In 
order to achieve more accuracy, Bartlett (1937) proposed multiplying LR by 
a constant, (1 + C/p)~ l , thus obtaining what is now known as the Bartlett- 
corrected test statistic: 

rp*_ LR 

where C is a constant of order n~ l chosen such that, under Ho, E(LR) = 
p{(l + C/p) + 0(n~ 3 / 2 )}. In regular problems, and under the null hypothesis 
LR* is \% distributed up to an error of order rT 2 ] see Barndorff-Nielsen and 
Hall (1988). A general expression for C in terms of log-likelihood cumulants 
(up to the fourth order) was obtained by Lawley (1956). 

One of our goals is to obtain the Bartlett correction term C for the test of 
7io : if> = against TCi : if> ^ if>^ for mixed linear models. This is done 
in Appendix A using Lawley's results; see (A.l). For simplicity, here we only 
give the expression for C when the ip^ = 0, which is common in practical 
applications: 

C = tr (V 1 [- l -M + \p - \ (7 + v) r T }) , (7) 

where tr(-) is the trace operator. Here, D, M and P are (m + 1) x (m + 1) 
matrices given by 1 

D = {(l/2)tr(E»"Efc)}, 
M = {t?{{X'^X' p )-\X?ti k X' p + 2X?±?X' p ))}, 

P = {t.{{X' p ^X' p ){X'^^ 

and r, 7 and v are (m + l)-vectors whose jth elements are tr((Xp T S _1 Xp) _1 
(X; T S^)), \x{D~ x A&) andtr((X T S- 1 X)- 1 (X T S^'X)), respectively. In our 
notation, is an (m + 1) x (m + 1) matrix given by 

A® = {(1/2) ti(±% k ) - (1/2) tr(E%) - (1/2) tr(£ j £ lk )}. 

Also, ±j = dH/duj, ± j = dH^/dujj = -E^E^E" 1 , E ife = <9 2 E/<%&j fe , 
fjk = dPH-i/dwjduk = -2S fc S j S- 1 - E^E^E" 1 and Xj = dX'Jdwj = 

-x^^xy^t^Xp. 



1 Note that we give the (j, k) element of each matrix. 
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It is noteworthy that (7) generalizes the result in Zucker et al. (2000, eq. (3)); 
their expression is only valid when the parameter under test is scalar and the 
covariance matrix for the random effects has a linear structure (and so does 
£, i.e., S = J^^jQj, where Qj are known matrices). Note that when £ has a 
linear structure we have £ j = Qj,Vj, = 0, Vj, k, and equation (7) becomes 

C = tr ( D - 1 {-iM + ip-i^}). (8) 

Additionally, when ip is scalar, our expression (8) reduces to equation (3) in 
Zucker et al. (2000). 

Cox and Reid (1987) proposed an adjustment to the profile likelihood function 
which can be used when the nuisance and interest parameters are orthogonal. 
The Cox-Reid adjusted profile log-likelihood function is given by 

iM) = W) -iiog{| -t H (£(V0) | } , 

where l^, is the matrix of second derivatives of £ with respect to 0. The 
corresponding likelihood ratio test statistic is 

LR CR {^) = 2 {VW-V(^ (0) )}. 

where ip is the maximizer of £ pa {ip)- 

The Cox-Reid test statistic is \% distributed under Ho up to an error of order 
rT x (like the standard likelihood ratio test statistic). DiCiccio and Stern (1994) 
defined a Bartlett correction to this test statistic which reduces the order of 
the approximation error to 0(n~ 2 ). The corrected test statistic is 

r d* _ LR CR 

LKcR - l + C*/p> 

where C* is a constant of order rT x such that, under Hq, E(LR C r) = p{(l + 
C*/p) + 0(n~ 3//2 )}. A general expression for C* can be found in DiCiccio and 
Stern (1994, eq.(25)). In Appendix B, we obtain C* for the test of TLq in mixed 
linear models; see (B.l). Here, we give the expression for C* for the case where 
^(°) = 0: 

C* =tr (V^-M + ip + 7 *r T }) , (9) 

where D, M, P and r were given above and the jth element of the vector 7* 
is tr(D _1 C < -^), with being an (m + 1) x (m + 1) matrix given by 

= {-tr^SjE" 1 !:,) + (1/2) tr(E^) + (1/2) tr(S fe S j7 )}. 

Our expression for C* generalizes the result in Zucker et al. (2000, eq. (4)), 
since their formula is only valid when p — 1. We notice that their formula 
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remains valid when the covariance matrix for the random effects has a nonlin- 
ear structure. As expected, (9) reduces to equation (4) of Zucker et al. (2000) 
when p — 1. 

The expressions we give for C and C* in (7) and (9), respectively, only involve 
simple operations on vectors and matrices. Therefore, they can be easily com- 
puted with the help of a programming language or software which can perform 
such operations, e.g. Ox (Doornik, 2006) and R (Ihaka and Gentleman, 1996). 
We note that C and C* only depend on X, on the inverse covariance matrix 
on the covariance matrix £ and its first two derivatives with respect to 

UJ. 



4 Numerical evidence 



In this section we shall present the results of a Monte Carlo simulation in which 
we evaluate the finite sample performances of the likelihood ratio test (LR), 
its Bartlett-corrected version (LR*), the adjusted profile likelihood ratio test 
(LR C r) and its counterpart in which the test statistic is Bartlett-corrected 

( lr cr)- 

The simulations were based on the following mixed linear model: 

Yij = ft + PiUj + (3 2 x 2 i + ftz 3 i + b 0i + buUj + e^, 

for j = 1, 2, . . . , Tj with Tj G {2, 3, 4, 5, 6, 7, 8, 9} and i — 1, 2, . . . , N. The values 
of tij were obtained as random draws from the standard uniform distribution 
U(0, 1); X2i and x 3 i are dummy variables. The fixed effects parameters are 
ft, ft, ft, ft- Also, b, = (boi buY ~ A/" 2 (0, G) with 



G 



U> 2 UJ 3 



(10) 



Additionally, the e^'s are independent from the foj's, and ~ A/" Ti (0, cj 4 / t J. 
We test Ho : ip = against Hi : ip ^ 0, where ip = ((3 2 ft) T . 

All simulations were performed using the Ox matrix programming language 
(Doornik, 2006). The number of Monte Carlo replications was 5,000 and the 
sample sizes considered were iV = 12, 24 and 36. The parameter values are 
p = 0, ft = 0.2, ft = 0, ft = 0, U! = 1, u 2 = and 0.25, u 3 = 0.5 and 1, and 
u 4 = 0.05. All tests were carried out at the following nominal levels: a = 5% 
and a = 10%. 

The null rejection rates of the four tests under evaluation are displayed in 
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Table 1. We note that the likelihood ratio test is liberal. For instance, when 
002 = 0, = 0.50, N = 12 and a = 10%, its rejection rate exceeds 20%. It is 
noteworthy that the three alternative tests outperform the standard likelihood 
ratio test. For N = 12 and N = 24, the two best performing tests are LR C r 
and LR* CR ] LR* is slightly oversized. For example, when uj 2 = 0, cj 3 = 0.50, 
N = 12 and a = 5%, the null rejection rates of LRcr, LR* cr and LR* 
are, respectively, 4.5%, 5.3% and 7.6% {LR: 13.0%). It is not possible to 
single out a global winner between LR C r and LR* CR . When N = 36, the 
Cox-Reid and the two Bartlett-corrected tests still outperform LR; here, LR* 
slightly outperforms the other two alternative tests, LR* CR being the second 
best performing test. 

Table 1 

Null rejection rates of the tests of Hq : tp = 0; entries are percentages. 











a 


= 5% 






a 


= 10% 




N 






LR 


LR* 


LRcr 


LR *CR 


LR 


LR* 


LRcr 


LR* C R 


12 





0.50 


13.0 


7.6 


4.5 


5.3 


20.8 


13.1 


9.2 


10.2 







1 


13.4 


7.8 


4.8 


5.9 


21.7 


13.5 


9.6 


10.8 




0.25 


0.50 


11.2 


6.0 


3.4 


4.1 


19.0 


11.2 


7.5 


8.5 




0.25 


1 


13.8 


7.9 


5.1 


5.8 


21.9 


13.9 


9.6 


10.7 


24 





0.50 


8.3 


5.6 


4.7 


5.0 


14.6 


10.9 


9.5 


10.0 







1 


8.5 


5.8 


4.9 


5.1 


14.6 


11.1 


10.1 


10.5 




0.25 


0.50 


8.6 


5.7 


4.8 


5.1 


14.8 


11.1 


9.6 


10.2 




0.25 


1 


8.7 


6.0 


4.8 


5.1 


15.0 


11.4 


10.1 


10.6 


36 





0.50 


6.4 


4.6 


4.2 


4.4 


12.8 


10.1 


9.5 


9.8 







1 


6.1 


4.9 


4.4 


4.7 


12.6 


9.8 


9.0 


9.4 




0.25 


0.50 


6.7 


4.8 


4.3 


4.6 


12.4 


10.0 


9.3 


9.6 




0.25 


1 


6.4 


4.7 


4.3 


4.4 


12.6 


9.8 


9.1 


9.4 



Figure 1 plots the relative quantile discrepancies against the asymptotic quan- 
tiles (N = 12, the smallest sample size, where the corrections are mostly 
needed). Relative quantile discrepancies are defined as differences between ex- 
act (estimated by simulation) and asymptotic (x|) quantiles divided by the 
latter. The closer to zero these discrepancies, the better the approximation 
used in the test. We note that the test statistics with the smallest relative 
quantile discrepancies are LR C r and LR* CR . We also note that quantiles of 
LR are approximately 50% larger than the respective asymptotic (x|) quan- 
tiles. 



5 An application 

We shall now present an application that uses a real data set. The data consists 
of a (randomly selected) subset of the data used by Crepeau et al. (1985). Heart 
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Fig. 1. Relative quantile discrepancies plot: N = 12, L02 = and W3 = 0.50. 

attacks were induced in rats exposed to four different low concentrations of 
halothane (group 1: 0% (control), group 2: 0.25%, group 3: 0.50% and group 
4: 1.0%); our sample consists of 23 rats. The blood pressure of each rat (in mm 
Hg) is recorded over different points in time (from 1 to 9 recordings) after the 
induced heart attack. The main goal is to investigate the effect of halothane 
on the blood pressure. 

At the outset, we consider a model where blood pressure grows linearly with 
time, possibly with different intercepts and slopes for each concentration of 
halothane, and with intecept and slope random effects to account for animal- 
to-animal variation. The mixed linear model can be stated as follows: 

Uij = A) + PiUj + 7o2G 2 i + 7o3C 3 i + 704^ + li2G 2 iUj + ^izGziUj 
+ 'juG&tij + boi + butij + 6ij, 

with i — 1,2, ... ,23 and j = 1, 2, . . . , r iy where yij is the blood pressure of the 
ith rat at time j, tij is the jth point in time (in minutes) in which the ith rat 
blood pressure was recorded, and G21 is a dummy variable that equals 1 if the 
ith rat belongs to group 2 and otherwise. Also, G 3i and G 4i equal 1 for groups 
3 and 4, respectively. We assume that bi = (&oi &ii) T '~ ^(0, G), where G 
is given in (10). Additionally, Af(0, (J4), the e^s being independent of 

the Vs. 

The maximum likelihood estimates of the fixed effects parameters are f3o = 
104.360, /3i = 0.004, 702 = -0.719, 703 = 0.203, 7 04 = -15.211, 712 = 0.022, 
7i 3 = 0.109 and 714 = —0.019. We wish to make inference on 7 12 , 713 and 
714. More specifically, we wish to test H : ip = against Hi : ip 7^ 0, where 
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■0 = (712, 7i3) 7i4) T - Note that under the null hypothesis, the mean slopes are 
equal for the different halothane concentrations. The adjusted profile maxi- 
mum likelihood estimates of 712, 713 and 714 are 712 = 0.020, 713 = 0.101 
and 7i4 = —0.030, respectively. The test statistics assume the following val- 
ues: LR = 6.522 (p- value: 0.089), LR* = 5.678 (p- value: 0.128), LR a = 5.287 
[p- value: 0.152) and LR* CR = 6.168 (p- value: 0.104). The standard likelihood 
ratio test rejects the null hypothesis at the 10% nominal level, i.e., it suggests 
that there are differences in mean slopes for different dosages. The three mod- 
ified tests (the Cox-Reid test and the two Bartlett-corrected tests), however, 
suggest otherwise, i.e., the null hypothesis is not rejected by these tests at the 
same nominal level. 

We now consider the following reduced model: 

Vij = A) + PiUj + 702^21 + 7o3C 3 i + 7o4G 4 i + b i + buUj + e i3 -, 

with i = 1, 2, . . . , 23 and j = 1,2,..., 7*. We now want to test H* :■?/>* = 
against H\ : ip* 7^ 0, where ip* = (702, 7o3 ? 7o4) T - Note that we are test- 
ing whether the mean blood pressures are equal across the different dosages. 
The fixed effects maximum likelihood estimates are (3q = 99.531, f3\ = 0.006, 

702 = —0.525, 703 = 2.318 and 704 = —13.357. The adjusted profile maxi- 
mum likelihood estimates of 702, 703 and 7 04 are, respectively, 702 = —0.823, 

703 = 2.079 and 7 04 = -12.573. We now have LR = 6.143 (p- value: 0.105), 
LR* = 5.174 (p-value: 0.159), LR a = 4.002 (p-value: 0.261) and LR* CR = 4.167 
(p- value: 0.244). All tests yield the same inference, namely: the null hypothe- 
sis is not rejected at the 10% nominal level, which agrees with the findings of 
Crepeau et al. (1985). 



6 Concluding remarks 

We addressed the issue of performing likelihood-based testing inference on 
the fixed effects parameters of mixed linear models when the sample contains 
a small number of observations. The standard likelihood ratio test is liberal 
(oversized), as evidenced by our Monte Carlo results. We obtained three al- 
ternative tests, namely: an adjusted profile likelihood ratio test, its Bartlett- 
corrected version and also the Bartlett-corrected likelihood ratio test. Our 
results generalize those in Zucker et al. (2000) in two directions. First, we 
allow practitioners to test joint restrictions on one or more fixed effects pa- 
rameters, whereas their results only hold for tests on a parameter at a time. 
Second, unlike Zucker et al. (2000), we do not assume that the covariance 
matrix of the random effects is linear when deriving the Bartlett correction to 
the profile likelihood ratio test. Our main results are stated through closed- 
form formulas that only involve simple operations on vectors and matrices, 
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and hence they can be easily implemented in matrix programming languages 
and statistical software. The numerical evidence we report clearly show that 
the proposed tests outperform the standard likelihood ratio test, especially 
when the sample size is small. It shows that the three alternative tests yield 
reliable inferences even for unbalanced data. In particular, the adjusted profile 
likelihood ratio test and its Bartlett-corrected version dramatically improve 
the type I error rate, especially when the number of observations is small. 
We strongly recommend that practitioners who wish to make inferences on 
the fixed effects parameters that index the mixed linear model on the basis 
of a small number of observations use the more reliable tests proposed in this 
paper. 
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Appendix A. Derivation of C 

Lawley's (1956) formula for C is 



where the indices r, s, t, u, v, w refer to the components of d = (ip , £ T , cj t ) t . Here, 
£ u> denotes summation over all possible combinations of the n+m+1 parameters 
in and w denotes summation over the combinations of the n — p + m + 1 
parameters in (£ T ,a; T ) T . We use the following tensor notation for log-likelihood 
cumulants: 



where C\ 





and 




{} r being the rth element of The notation used for derivatives of cumulants is as 
follows: 

(In what follows, we shall use similar notation for derivatives of matrices formed 
out of cumulants.) Note that —K rs is the (r, s) element of Fisher's information 
matrix; the (r, s) element of its inverse is denoted by —k ts . We use indices a,b,c,d 
in reference to the components of ip, indices /, g for the components of £, and indices 
j, k,l,o for the elements of u. Further notation used here is given in Sections 2 and 
3. 

The first-order derivatives of the log-likelihood function in (6) are 



0£(0;Y) 



— itr(S -1 Sj) - -z T £ J z + V T ^' T S- 1 z. 



^ 2 v J/ 2 

The second-order derivatives are 

j k 

+ X£ T 5^ + Xj T S fe )z - iz T S^z, 

where 

Jr' = _^ = 2X(x T E- 1 x)- 1 x T s fe x(x T s^ 1 x)~ 1 x T i; i x' 

^ duj k ' y ' p 

- X(X T 'S- 1 X)- 1 X T, E jk X' p . 

Additionally, the third-order derivatives are 

d 3 m;Y) _ ~ T . ■ ~ dH($;Y) _ ~ /T • • ~, 

g^Y) g 3 l(l;Y) WlZl = x T ^ fc fY-Xn 

dipdcojdujk ]k k 3 p ' 

; ^ Y ' J( tr(S Zfc S,) + tr(S fe E^) + tr(S'S jfe ) + tr(E _1 Sj M ) + z T E jfc ,z) 



dujdujkduji 2 
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where itjki = d'Tijk/duoi and X'- kl = dX'- k /du)i. Finally, the fourth-order derivatives 
can be shown to be 



d 4 £{tf;Y) 

d 4 £(tf;Y) dH($;Y) dH{$;Y 



2X'J± j xL -x'J£ jk xl, 

0. 



dt/jdtp T d^fd^g di/jdip T d^fdujj dtpdip 1 <9£/ 
Taking expected values of second, third and fourth derivatives, we obtain 





= E 




= E 




= E 



- x y x 





= -X T S~ 1 X, 








= 0, K£ u . Uk = 



In similar fashion, it follows that 

= X T ± j X' p i/>, K^. = 0, 

X S J ' fc Xp ^, KiW£ f £ g = K^ fUlj = 0. 
Additionally, 

/e yfc = -2 ^(S'SfeS- 1 ^) + X - tr(SJ'S, fc ) + i tr(S fc E y ) + i tr(E'i) jfc ). 

Consider the following matrices formed out of minus Fisher's information inverse: 
K** = K~l, = (K^ - Kj u K£Kz u )-\ = K' 1 + K^ 1 K iu IC^ Kj u K^ 

and = K^ T = -K^ 1 K iu K uuT , where the jth column of is K^ u . and the 
(j, fc)th element of K ww is Kj k - It can be shown that 

( K W>)j = -Xp&X'p, ( K 4>ip) jk = -zx'k'EPx'p - x p T j: jk x' p , 
(K^) j = -X T VX, (K (u .) k = X T & k X' p $ + X T ± j X' k iP, 

{Kji) k = -trtS'Sj-E^Sfc) + -ti{± j ± lk ) + -ti{±% k ). 

It follows from the orthogonality between tp and (£ T ,o; T ) T that K a ? = K a i = 
( K af)jb = (naj)fb = 0. Also, Kjf a = Kjf ab = 0. Hence, 

C\ = (labcd + labfg + Lbfj + Lbjf + Lbjk + If gab + Ijkab) > 

where ^ ranges over all parameter combinations induced by the indices a, b, c, d, f, g, 
j, k. It is possible to show that l abcd = l abfg = l abfj = l abjf = l fgab = 0. Thus, 

Cl = ^ ] (Jabjk + Ijkab) = ^ ] ^ ^^^abj'fc ~~ i^abj)^ + K ^jfcob^ • 
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Since K abjk = (K a bj) k = Kjkab, Ci reduces to 
As for C2, we have that 

C*2 ^ ^ (Jabcdjk ~l~ labjkcd 4" Ijkloab 4" Ijkablo 4~ Ijkabcd 4~ Ifjgkab 4~ ^fjkgab 4~ Ifjklab 

4~ Ijkfgab 4" ^jkflab 4~ ^jklfab 4~ ^jfabgk 4~ ^jfabkg 4~ ^jfabkl 4~ Ijkabfg 4~ Ijkabfl 4" ^jkablf) 

E( ^- ab cd ^jk , ^ ^.ab fg jk 

^ — ft, ft ft ruabj K cdk "T —ft ft ft ^-abj^fgk 

- K ab K fj K M K abj (^2(ft /fc )z - |k /w ^ + ^^^W^^ifo - 2(K W ) )}. 

Therefore, C reduces to 

C = ^ ] I — —ft ft^ K a bjk + ^ K K ^ KabjKcdk ~ ~K. K? K K a },j (ft; fc — 2(ft/ Q )fc) 
+ K ab K fk K jl K abj (2(Kf k )l - |k /W ) ~ ^K ab K f9 K jk K abj K fgk y 

We now arrive at the matrix expression given by 

C = tr (V- |-lM + ip - Q p - 5 + \ r T }) . (A.l) 

Here, p, 5 and 77 are (m+l)-vectors whose jth elements are, respectively, tr^^A^)), 
t r (K^ T j g(i)) an d tr(-K^(X T t,iX)). In our notation, fiO') i s a matrix that con- 
tains the m + 1 column vectors (1/2 X T 1t? k X' p + 2 X T i^X;(.)Y> and A^ is defined 
in Section 3. For the test of Ho 4> = 0, C reduces to equation (7). 

Appendix B. Derivation of C* 

We shall now obtain C*, which is used to Bartlett-correct the adjusted profile like- 
lihood ratio test statistic. DiCiccio and Stern (1994, eq. (25)) give the following 
general expression: 

C* = \ \r ru r st K rstu - ft"V st (ftv st ) u + (n ru K st - v "V st ) ( Krs ) tu 

_ (l K ru T st T vw + }_ K ru T s Wr tv _ ^ru^w^ + (^ru^t^w + ^u^sw^tv 

- u ru K sw v tv ) K rst (k uv ) w - (n ru K st K vw - v "V st u ™) {K rs ) t {n uv ) w 

- [K rU K SW K tV - V "V SW U tV ) (K rs ) t (K uv ) w 1 , 
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where v rs = n rs — t ts , r rs = n n sa a ab , c ab being the (a, b) element of the inverse 
of . From the orthogonality between tp and <fi we have that t? 9 = r jfe = T*i = 

T af = T aj = q Alg0) T ab = ^ab _ 

C* = |— K ad K b ° Kafrcd — K rU K ab (K rab ) u + (k TU K St — V TU V S *) {n r s)tu 

~ (4^ K ^ + 2^ ^ K ) ^rab^ucd ~\~ K K K K rab {^uv) w 
+ ( K ru K tv_ u ru u t^ K s WKrst j. 

We have that K ru K to - z^"V to = K "V to + K to r rM - r"V to and (K 6d ) fc = K 6dfc . Hence, 

^ ^ (n K — V V ) K K rs t (k uv ) w = ^ ^ K K K? t^ajc^bdk- 

Since n abcd = K abc = n fab = (n ac ) bu = (K a f)tu = (n aj ) tu = 0, it follows that C* 
reduces to 

C = ^ ^ ^ K Kabjk + K K abj K cdk + K K? K ^abj{ K kl)o 

-\- K K 9 K ab j Kfgfc -\- 2k K K ab j ( K fk ) [ j- . 

We then arrive at the matrix expression 

C* =tr j-M + ^P +(/>* + 2 <T)r T j^ + r T i^r/*, (B.l) 

where the j'th elements of the vectors p* and 5* are, respectively, tr(if a "*'C^^) and 
t r (K^ T F^), and the /th element of the vector n* is tr(K^G^). Also, C^') is 
defined in Section 3, is a matrix that contains the m + 1 column vectors 
(x T E^X£ + X T S^) ^ and is the (n - p) x (m + 1) matrix whose j'th 

column is the /th column of — X T ~E^X. For the test of Hq : = 0, C* reduces to 
equation (8). 
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